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ABSTRACT 


A semi-empirical method for modelling the loss of electron fluxes in the earth’s mag- 
netosphere was developed. An equation for the integral-energy omnidirectional electron flux 
as a function of time and magnetic field strength was derived from pitch-angle diffusion 
theory. 


This flux equation was the basis for a computer data-fitting program written at the 
Air Force Weapons Laboratory (AF WL) to fit the AF WL Trapped Electron Daia Base. The 
program utilized a Jeast-squares fil and incorporated random variations of the characteristic 
exponential loss times about their initial values. An improved table of initial lcss times was 
compiled for use with the program. 


The derived flux model showed substantial agreement with the empirical data base. 
Representative plots of computed flux over raw data are shown for L-values of 2.3 and 2.4. 





1. INTRODUCTION 


Objective 


Modern technology, both civilian and military, depends greatly on the use of satellite 
systems. It is thus important to the Department of Defense to understand and predict the 
environment in which satellites may operate, in order to determine survivability and vul- 
nerability requirements. One part of that environment is the electron fluxes which may exist 
in the trapping regions of the earth’s magnetosphere. Such fluxes may arise from natural 
sources or from injection by high-altitude nuclear bursts. Because such fluxes have great 
capability for inflicting damage to sensitive satellite components, even at great distances 
from a burst, understanding their behavior is vitally important. 


There are a large number of theoretical treatments in the literature which describe 
various electron loss mechanisms from the earth’s magnetosphere. However, the empirical 
calculations of actual electron losses are often conflicting, and frequently are restricted in 
their areas of coverage. 


Since the cessation of high altitude nuclear testing in 1962, loss studies have been limited 
to electrons injected into the magnetosphere by natural magnetic disturbances such as solar 
storms. Thus, much of the theory which has been developed since 1962 has had restricted 
opportunity for empirical testing. 


The objective of this study is to produce an improved algorithm to calculate losses 
of electrons from the magnetosphere. The algorithm is semi-empirical, because it utilizes 
the Air Force Weapons Laboratory (AF WL) Trapped Flectron Data Base (Ref. 37) measure- 
ment of electron fluxes (following the old high altitude nuclear bursts) to determine the 
characteristic loss times of those electrons. 


The loss times thus determined and the algorithm for ux “decay” will then be incor- 
porated by AF WL into the SPECTER computer code (Ref. 7-9) for improved flux calcula- 
tions of satellite environments. 


Scope 


This study is limited to investigation of one loss mechanism, pitch-angle diffusion, using 
a theoretical formalism developed principally by Shulz (Ref. 47,48). The study has also been 
necessarily limited to those regions of space and those energies covered by the satellites from 
which the data base was complied. The loss equations thus developed in this study are only 
implicitly functions of energy, E, and L-coordinate (explained in Chapter II). 


The sheer volume of data and the limited time for this study also limited the testing 
of the algorithm developed to a representative subset of the entire data base. 


Approach 


A specific solution of the pitch-angle diffusion equation is derived using the method 
proposed by Shulz (Ref. 47). This solution is then fitted to the electron flux data using 
a least-squares fitting program written by AFWL, with the coordination of this author. 
The solutions are performed for integral fluxes (above threshold energy of the satellite 
detector) and for constant L-value. The solutions are thus explicitly functions of magnetic 
field strength, B, but implicitly functions of / and L. 


Use of the AFWL/NTCTS computer facilities was essential to this study. Well over 72 
hours of actual computer- processing time were used, and over 1000 computer-produced flux 
plots were examined. 


Presentation 


A brief Background section is presented, which covers some physical relations and 
terminology used in the study of space physics. This section is used to establish uniform 
symbol notation, and to provide a common starting point for the reader unfamiliar with the 


specialized language of space physics. It may be skipped by the reader more familiar with 
the topic. 


Chapter III outlines the theory used to develop the flux equations which were fitted to 
the data. 


Following the theory chapter is a chapter outlining the method of data analysis. 
Descriptions are given of the data base and of the AF WI computer codes used in this study. 


Representative results of the study are presented in Chapter Y, in both tabnlar and 
graphic form. Again, the sheer volurne of data makes inclusion of all plots and fitted values 
prohibitive. 


Conclusions and recommendations for improvements for further study are presented in 
Chapter VI. 


The Bibliography includes 65 sources and provides 2 comprehensive summary of existing 


literature of relevant topics about radiation trapped in the magnetosphere and of pitch-angle 
diffusion. 
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H. BACKGROUND 


The equation cf motion of a charged particle in magnetic, electric, and external gravita- 
tional fields is: 


d*7 a dr = = 
m de = mg+a(o x it) + a8 
where 

Y = Fibro, Uo): a — 
m = particle mass, 
gj = gravitational constant, 
q = charge, 
B= magnetic field, 
iE = total clectriec field (Ref. 11:23-24.) 


In a uniform magnetic field, a charged particle will move in helical fashion along and 
around the field lines with a cyclotron radius, p,, (also called gyroradius or Larmor radius) 
about its “guiding center”. The particle’s motion may thus be separated into its rotation 
around the field line and the motion of its center of rotation, or guiding center, along the 
field line (Ref. 11:24-25). 


The period of cyclotron rotation, 7, is defined by (Ref. 39:5): 


Fa a 


where m is the relativistic mass: 


m = m7 = my| ————_— 
1 — (2 /c?) 


and where v; is component of velocity perpendicular to B. The cyclotron frequency is 
simply 2x over the period (Ref. 39:5): 


of qB 
ie 


ve S = 


The angle between a particle's local velocity, v, and the magnetic ficld, B, is its pitch- 
angle, a, which is defined as: 


v 


I fv 
a = arceos(—) = aresin(“), 
Uv Uv 


In a uniform static magnetic field, vy is constant and vy is constant, and hence a is constant. 
The particle will then move with a uniform circular motion around the field line and uniform 
rectilinear motion along the field line, resulting in a helical motion about the line, sinec 


dp/dt = qovxB = G0. 
and 
mi, = gv X B), =o), 


a) 


dy =~ 0, 
uy = conatant, 
and 7 
ma; = qx B),, 
a, = (qvyB)/m = _ constant, which is 4 constant centripetal arceleration 


(Ref.39:4-7; Ref. 11:23-36). 


The concept of the first adiabatic invariant arises naturally from the guiding center 
approximation. The first adtabatie invariant, M, (also called the relativisi:ic magnetic 
moment) is defined by: 

Pi 


2mo 





M= = conatant 

where p, is the particle's perpendicular momentum in the guiding center ap oroximation. 
The assumptions implicit in calling M invariant are that the spatial variation >f /3 is sinall 
compared to p- and the time variation of B is small compared to r, (Ref. 39:1 1-23). 


If a particle has a constant velocity along a field line in the guiding certer frame of 
reference, then 
sin- a(a) _ sin’ a; = 
“Bae = ae = constant 
where a is the arc length along the field line and z is any point on the field line. ‘This assumes 
that the particle’s kinetic energy remains constant as it follows the field line. If a particle 
enters the field at point B; with pitch angle a,, its velocity along the field line is 


B ee 
vy (4) = v cosa(s) =v @ \/ 1— a Sin” Qj. 


If the field is increasing in the direction of vy, then the “mirror point” of the particle, where 
its parallel velocity is zero, is 
B, Bla 
Be = aa = se 
sin-a; sin’ a(s) 





and the particle has a local pitch angle of 90° at that mirror point (Ref. 39:34 -42). 


When the magnetic field has a geomctry like that of a dipole field, increasing at the 
poles and decreasing in the midpoint, then the particle is in the so-called “magu:etic bottle”, 
and is trapped between the mirror points. There is some minimum B value between the 
mirror points which is called Bo, the equatorial field strength. Using the giiding center 
approximation and approximating the earth’s magnetic field by a dipole fiel!, a particle 
trapped in the earth’s field has three distinct motions. It circles rapidly about :. field line, it 
bounces along the line between mirror points, and it drifts slowly in longitud: around the 


earth. All three motions take place with different speeds, so they are distinct (Ref. 11:25; 
Ref. 29:34-65). 


From the mirror point definition, and from the velocity equations, we sce that 


Py 
y= —-— 
2 mo Bm 
and that 
B B. 
vy) (s) = sina; = = i 


Note that the particle’s mirror point depends only on its “injection” point and its 
injection pitch angle, not on its initial energy or velocity. A mirror point is a consequence 
of the field alone, and all particles injected at the same point with the same pitch angle will 
mirror at the same point on the field line. Of course, this is only urue if no external forces 
are acting. 


The parallel velocity, vq, of the particle will be the maximun at Bo, wher> Ba) is the 
minimum. Thus, a trapped particle spends most of its time near the mirror po-nts, and the 
least amount of its time transiting the equatorial field regions (Ref. 29:34-44). 


The bounce period, n, of a trapped particle is generally much greater than its cyclotron 
period, r,, and is defined by: 





aoa [~ ds =?/™ da 
Se coe (8 (6) ) 
where 8m, and 4m, are the mirror points on the field line. 


As a particle in 1 dipole field like the earth’s bounces along the field line, it also drifts 
perpendicularly to the field linc, due to external forces, field gradients, and fie'd curvature, 
as well as other effects such as time-dependent field changes. This drift is slow compared 
to the bounce period. As the particle bounces and drifts, it traces out a surface between its 
mirror points and around the earth, called a drift shell (Ref. 39:9-19). 


The concepts of the second and third adiabatic invariants arise from the above behavior. 
If the forces acting on the particle remain almost constant over its bounce period, the second 
adiabatic tnvariant, J, is defined by 


sap mdemal py ds 


where py is the momentum component parallel to B and ds measures arc lenyth from the 
equator (Ref. 45:12; Ref. 6:3-31). If the forces acting on a particle remain almost constant 
over its drift period, the third adtabatic tnvartant, &, which is the magnetic flux enclosed by 
a drift shell, is defined by 


b= fie e di 


where Ag is the magnetic vector potential and where the integration is perfcrmed over 2 
curve, 8, which lies in the drift shell. The third invariant, ®, is defined and c »mputed for 
the drift shell of the guiding center, with constant field, and not for a drift shell which a 
particle may physically trace out under short-term conditions (Ref. 45:12; Ref 39:76-79). 


The guiding drift shell of a particle may also be referred to as an “invari:int surface,” 
which is com posed of field lines which end at the mirror points. The three adiaba: ic invariants 
uniquely define an invariant surface (Ref. 6:3-31 to 3-32). 


The more common set of parameters used to define a particle’s position in the mag- 
netosphere (or to define an invariant surface) is the B-L coordinate system. The 8 parameter 
is magnetic field strength, and L is the Mcllwain L parameter defined by L =:o9//tg¢ where 
Rj: is the earth’s radius (= 6371 km) (Ref. 6:3-33) and (Ref. 39:53). 


The value of L, therefore, is equal to the distance, in earth radii, of t}e equatorial 
point on a field line in a dipole field. If the ficld is not symmetric (not a pe: fect dipole), 
the invariant surface is not so well defined. llowever, for most field lincs in -he Lrapping 


a | 


region of the earth’s magnetosphere, L varies by less than one percent along the line, so 
the B-L system is adequate. The less-than-exact symmetry of the earth’s dipole field and 
outside forces may cause particle drift to deviate from a perfect azimuthal course. This 
means that its / value will vary in the course of a drift period. This variatior of L values 
is not significant below L = 3, but may be so above that value. The average j-value of all 
intersecting invariant surfaces is called an L-shell. This L-value defines a set of surfaces along 
field lines, which may end at differing mirror points. The /-shells are considerec. to intersect 
the earth’s surface, even though mirror points do not extend through the atmo: phere to the 
surface (Ref. 6:3-33 to 3-35; Ref. 39:53). 


The various mechanisms which operate to cause trapped particles to chanze /-shell or 
to be lost from the trapping region will necessarily violate one or more of the adiabatic 
invariants. For example, pitch-angle diffusion violates M or J, or both; and racial diffusion 
violates © (Ref. 49:48). These concepts are discussed further in Chapter III. 


IM eels rete © Tay, 


While the motion of particles trapped in a magnetic field is well understood, the 
mechanisms of losses from magnetic fields such as the earth’s magnetosphere are Icss well 
understood. It is generally agreed, however, that pitch-angle diffusion into the loss cone is 
one of the predominant mechanisms for removal of charged particles for mid-range /-values. 
The loss cone angle is the lower limit of pitch angle for trapped particles. Any farticles with 
smaller pitch angles will mirror in the sensible atmosphere and will be lost by atmospheric 
scattering. 


Roberts (Ref. 41) has noted that a pitch-angle scattering mechanism must be extant 
for pitch-angle diffusion to occur. Such a mechanism would necessarily violate one or more 
of the adiabatic invariants. 


It is not the purpose of this study to ascertain the true source mechani3m of pitch- 
angle diffusion. Such a study is beyond the scope of this paper. Indeed, the search for 
this mechanism has been going on for decades. Roberts (Ref. 41) has postulated that the 
“perturbation forces” causing such diffusion may result from turbulent ambient electric or 
magnetic fields or from collisions with other trapped or non-trapped particles. Lyons (Refs. 
2-28) has formulated extensive theory based on resonant interactions of so-callid “whistler” 
VLF waves with the trapped-particles’ gyrofrequencies: “cyclotron-resonance”. Regardless 
of the actual physical mechanism behind pitch-angle diffusion, its general treatment is 
mathematically the same, and the physical results are the same. (Ref. 41:308) This study 
follows the methods of Roberts (Ref. 41), Shulz and Lanzerotti (Ref. 45) ano Shulz (Ref. 
A7) in assuming that pitch-angle diffusion is an operative process. One of the purposes of 
this study is the formulation of the equations necessary to validate that assumption against 
experimental data, and to perform that experimental validation. 


It should be emphasized that, while pitch-angle diffusion is assumed in this paper to 
be the dominant loss mechanism for L-values and altitudes considered, it is not the only 
loss mechanism. Several investigators have formulated radial diffusion (cros:-L) theorics 
(for example, Walt (Ref. 58) or Tomassian (Ref. 54)), coulombic or collisioa-scattering- 
diffusion theories (for example, Wentworth (Ref. 60)) and multiple diffusion theories, such as 
combinations of pitch-angle, energy, and/or radial diffusion (for example, Walt (Ref. 58) and 
Lyons (Ref. 26)), all of which show some agrcement with experimental data. Ia particular, 
atmosphcric scattering is obviously a dominant force at very low L values (Ref. 30). 


The introduction of multiple loss-emechanisms makes explicit solutions of :ny diffusion 
equation extremely difficult. To simplify the problem and render it amenable to the method 
of Shulz (Ref. 47), this study assumes that “as a rule, radial diffusion enables the radiation 
belts to become populated from an external source (or rearranges particles injected by an 
internal source), while pitch angle diffusion causes particle loss to an atmospheric sink” (Ref. 
45:48). 


Pitch-angle diffusion, while simple in concept, is complicated in detail. Roberts’ (Ref. 
41:307-337) treatment of the general mechanism is particularly descriptive: 


At the magnetic equator, a particle’s pitch angle, ao, is determined by 


P 


Jaren VU) 


(1) Z == cOSGy = 


“q 


where Py, and P,, are the parallel and perpendicluar components of the particle’s 
momentum with respect to the field line at the equatorial value of By. Roberts next defines 
2 
Pig 
omg 
Pye 
2m 





ilo = 


(2) 
WL = 


where mo is the rest mass. These are defined since the first adiabatic invariant, M, is 


proportional to W, |: 
W 
(3) M = ab 
| Bo 


x, 
and, if Wj), < W,,, the second adiabatic invariant, J, is approximately prc portional to 


Wi: 
rpW\, 





IR 

(4 ; 
where ; 

Tp = the particle’s bounce period, 


relativistic mass m 


y= —_ = 


rest mass mo 


A particle’s path may be defined by plotting its values of W1, and Wy, cs in Figure 1. 
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Fig. 1. Loss of a trapped particle by @ random walk into the pitch- 
anzle loss cone. (Ref. 41:306) 


Note that lines of constant equatorial pitch angle pass through the origin, because 





Any force which violates either M or J, the first or second adiabatic irvariant, will 
produce “diffusion” of the particle’s pitch angle as shown in Figure 1 as W;, and Wy, 
change with time. When a particle’s equatorial pitch angle diffuses to the vali e of the loss 
cone angle, 

Bo 
Byov 


where Boo is the value at 100 km altitude, then it is lost by atmospheric scattering. “The 
loss cone serves as a ‘sink’ for particles” undergoing pitch angle diffusion (Ref. 41:307). 





(5) 2, = cosa, = \/1- 


Roberts makes two very important points about pitch-angle diffusion as a oss-mechan- 
ism: 


(1) with no source, the entire radiation-belt would be depleted of particles, 


(2) the loss cone approximation is just that, since the atmosphere is not sh::rply defined 
at 100 km. However, if atmospheric pitch-angle scattering is not the primary focus, 
the loss come approximation may be useful. “Naturally, when the loss-:-one concept 
is used, detailed agreement between theory and experiment cannot be expected in 
the region near the edge of the loss cone” (Ref. 41:308). 


Both Roberts (Ref. 41) and Shulz and Lanzerotti (Ref. 45) develop the pitch angle 
diffusion equation from the Fokker-Planck equation using a particle distribution function, 


f. 


The Fokker-Planck equation is 1 formalism which arose from the study of Brownian 
motion, and which is used frequently in transport theory (Ref. 42:308). The characteristic 
Fokker-Planck equation for trapped particles (which ignores radial diffusion) i: 


of__ 1 O{ (ae) J, 1 0 a 
dt %wpdOE a(S ) + Oz [2 Ty) Des Oz 


where the first term, subscripted v, represents non-stochastic (mean) energy loss to the 
atmosphere, the second term represents pitch-angle diffusion, and the third ter n represents 
range-straggling (energy diffusion). The term T(y) is defined below, and the Der and D,, 
terms are the characteristic diffusion coefficients; 7 is relativistic mass ratio (Ref. 45:55-58). 





1 @ af 
eee ae Doe ee 
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Roberts uses a distribution of particles in a “tube” of force about a fir ld line, and 
Shulz and Lanzerotti use a phase-space density dtatributton function which ‘s essentially 
equivalent to Roberts’. Since this study follows the methods of Shulz (Ref. 4.°) and Shulz 
and Lanzerotti, a discussion of phase space is necessary (Ref. 45:15-22). 


Any moving particle can be described by specifying its three position coo-dinates and 
its three canonical momentum components. This completes a six-dimension “phase-space” 
in which a particle can travel in time. 


If there exists 1 system of a large number of particles in phase-space, thc system can 
be described by a six-dimensional distribution function (4, dy; t) where 


P; (¢ = 1,2,3) are canonical momentum components, 
Ti 


(: = 1,2,3) are position coordinates, and t is time. 


Thus, f d° Pd?q is the number of particles instantaneously occupying 6D volume 
d* Pd°3. According to Liouville’s theorem, the phase-space volume containing the system 
of particles moves incompressibly through phase-space. 


Since P is an awkward quantity to deal with physically, Shulz and Lanzerctti note that 
P transforms to the more familiar p as 


(8) P=p+ (t)A 
where A is the electromagnetic vector potential. 


Hence, S(F,4,t) = f(p,7,t) since ¢ = 7 and the P to p transformation has a unit 
Jacobian. The position-momentum distribution function f(p,7,t) defines the particles oc- 
cupying the 6-D volume d°pad°F at any given ¢. The pitch angle diffusion equaticn is given in 
terms of a phase-space distribution function which is numerically equal to the more easily- 
definable position-momentum distribution function (Ref. 45:15-22). At constant energy and 
L-shell, and under the action of some source S, the diffusion equation in phase-space can be 


written oe P a 
f l fy. s 
= —_ = — T Dez ——\ +585 
(") Ot z (ges t Ty) Des Ox 
where T(y) is the quarter-bounce integral path length function of y = 1-2? where 
= cosaq. The quarter bounce integral gives the length of the trajectory of « particle (in 
units of Ro, the distance to the equatorial crossing) from the equator to the rairror point. 
The exact definition of T is : 
: / ™ da 
T= -— 
9 cosa 


where 8 is the distance along the field line and a is the local pitch angle (Ref. 10:4029-4030). 


D,, is the bounce-averaged difusion coefficient as a function of z (Ref. 41 and Ref. 
45). The boundary conditions are that f = 0 at zr = z, (some cutoff value) ancl f(z = 0) is 
finite. 


Shulz and Lanzerotti (Ref. 45:162) rewrite the above as 


el) min 22 eer 
Zz E y ae T(y) Oz E 


and further approximate that the second term is negligible for 2 << 1 since r? + y? = 1. 
If T(y) ss T(1), the second term disappears and the equation is a diffusion equation in 
cylindrical coordinates. (These approximations are not used in this work. This jevelopment 
is used to aid in understanding the exact solutions which follow from Shulz .Ref. 47)). if 
D,2 and S are then independent of z, and if f(z, t) = X(z) T(t), the eigenfurctions of (8) 
are Bessel functions of order zero. The general solution to (8) would then be 





1 @ 
(g) Bais l|sp., 


9) H(2,t) = Jao(2) + Do an(t) Jo («. > 


-vhere f., is the steady-state solution 


eo rat)= (Te) (ZY 
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The x, are the zeroes of Jo (n = 1,2,...) and the a,,(t) vary as e7!'/'=) where rz = 
z? /(D*x?,). Thus, the pitch angle distribution (and hence the directional or omnidirectional 
flux) is shown to be the sum of a steady state and higher order eigenmodes (Ref. 45:160-168). 


The steady state can be thought of as the “normal” or “quiet-time” value of flux which 
exists in equilibrium with the source S. An “injection” of particles would theo result in a 
perturbed distribution function with several eigenmodes, each with a charactcristic decay- 
time t,. The higher order modes decay faster, and eventually only the fundainental mode 
would remain, which would decay exponentially to reach the steady state (after “infinite” 
time) (Ref. 41 and 45). There is ample evidence that such a process does indeed occur. 
Roberts (Ref. 41) cites Explorer XV data for the 28 October 1962 Russian explosion as an 
example. Rosen and Sanders (Ref. 43) also note that decay is faster immediately after solar 
magnetic storm activity than during quieter periods. The similarities in Figures 2 and 3 
illustrate the way in which a decay of higher eigenmodes (Figure 2) can approximate the 
temporal evaluation of an actual electron distribution (Figure 3). 


The primary difficulties with the above discussion arise from the approximations that: 
(1) Ty) = T(1) = constant 
(2) 2° <1 (ory 1) 
(3) Dz, varies little (or not at all) with z. 


Roberts (Ref. 41:311) uses the argument that the full spiral path (and hencc T(y)) varies 
only by a factor of 1.4 as z varies from 0 to 0.9 and by 1.9 as z varies from #) to 1. Shulz 
and Lanzerotti (Ref. 45:163) use the simplifying assumptions that T(y) = T(1) and D,, is 
independent of z. Both use the assumption that 2° < 1. 

The present study requires that z be allowed to vary from 0 to z, (up to near 1) in 
order to adequately examine the AF WL trapped electron data base (Ref. 37). In addition, 
if z is not small, the approximation that T(y) = T(1) is a poor one, since y is not close to 
|b 

The function T(y) can be shown, within 0.57% to be approximated as (Ref. 10:4030) 

(11) T(y) = T(0) — [T(0) - T(1)]y3/4 
where T(0) = 1+ (inte - v3))/(2v8) mz 1.3802 and T(1) = (xv2)/6 ss 0.7405 (Ref. 
45:19). 


Obviously, if z approaches 1, y approaches 0 and the T(y) = T(1) approximation 
becomes invalid. 


The three limiting approximations are removed in the treatment of Shulz (Ref. 47) by 
introduction of a new “canonical” variable z such that: 


z ] 
(12) z= ay) = fatty yas! = fy TY) ay 
; ¥ 
and a corresponding diffusion coefficient 


(13) J, == (z Ty) De. 
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Figure 2. Postulated Decay cf Higher Flux Eigenmodes to the Steady Strate by Pitch- 


Angle Diffusion (Ref. 41:313). 
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Figure 3. Temporal evolutien of 1.9 Mev Electrons from the Eussinn 3 Burst (Ref. 


41:319). 


Using (11) in (12) and performing the integration: 


is / y[ro- IT(0) — rca 


ae / row =(re— rly ay 


1 3 4 1/4 
~~ sT(O)y — SITO)-TUly I, 


royr—y2) — (70) - TAL 9/4 
2 . 11 
where Z(0) = 18/35 and 7(1) = 0 (Ref. 46:5213). 


(14) ne eS 


Using (7), (12) and (13), the diffusion equation becomes: 
_ aie 
(15) -~ = —'!D,,; ale oO. 


Sholz (Ref. 47:8) states that (15) is a “canonical” diffusion equation in that there is no 
Jacobian factor ‘which fails to commute with D,, as was the case in (7). 


The assumptions which now must be made are: 


(1) D,, is some “suitably simple function of z” (Ref. 47:6) so that ezact eizenfunctions 
gn(z) may be specified, 


(2) S, the distributed natural environment source, is independent of z. 


Assumption (1) is the limiting assumption, since if the function of z is very :omplex, the 
diffusion equation becomes extremely difficult to solve. Assumption (2) is both simple and 
reasonable, however. The source of the natural environment must be close tc constant in 
order to be the “driver” of quiet-time equilibrium. It must also be distributed ‘airly closely 
to the steady-state quiet time distribution for the same reason. 


Shulz (Ref. 47:7) states that even if D,, is not of a functional form tc yield exact 
eigenfunctions, there may be a D,, which resembles D,, “closely enough” aud for which 


exact eigenfunctions j,,(z) (resembling g,(z)) are known. 


The following derivation of the exact eigenfunctions and the correspondinz omnidirec- 
tional differential-energy flux arises directly from Shulz (Ref. 47). 


The first assumption by Shulz is that 


(16) Dx — (=) De .. 


where ¢ is some number, not aecessarily an integer, less than 2, and where ),- is the value 
of D., at some z_ < 16/35 where f vanishes. Roberts (1969) also uses this assumption. 


This form of D,, allows a basis set of orthogonal eigenfunctions to be s] own for the 
interval of interest: 0 < z < zy. Following tle notation above, we seek scme ),, and 


eigenfunctions g,(z) which “closely resemble” the true functions and which sa isfy 
(17) Da 2 Jnl) iol=)) + ind, (2) = 9 


for 


with the substitution 








where e 
€ =1- 9? 
l-o 
c= 9 
and 
\ale? 
i= = 2 
Dee (et 


then (17) becomes Bessel’s equation. The exact eigenfunctions of (17) are giver. by 
2-0 


(18) - = & = (: abe 
Gn (2) = Ys ee Ze Ja Kon ze 


where Ky, are the nth zeroes of the Bessel function J, of the order v, where 








The g,, given by (18) are normalized so that 


(18) a / In Im dz = bnm. 


The cigenvalues of (18) are given by 





(19) a & ye Pee Be 


Shulz and Boucher (Ref. 48:6) state that, since particles are not lost at z = 0 (but 
rather at z = z,}, the diffusion current must go to zero at z = 0. A diffusion current at 
z = 0 would imply diffusion of pitch angles (and hence mirror points) into the equator, 
which is not reasonable. The condition of no diffusion current at z = 0 co-responds to 
lim,—9 Dee qo) = 0. 


Shulz and Boucher (Ref. 48) use the series expansion of J, from Abramowic: and Stegun 
(Ref. 1) to show the limit as z approaches zero. This author shows the limit directly for the 
special case eigenfunction (J_1/2) which is used in this report. The derivatior is shown in 
Appendix A. 


In order to arrive at the differential-energy omnidirectional flux, Shulz (Ref. 47) assumes 
that f(z, t)—> f(z) as too; that is, the solution distribution function has some steady-state 
value. This also assumes that S is constant in time, but not in z or z. The f_, distribution 
is necessary in order to expand f as ascries of cigenfunctions in the flux equation. To arrive 
at an expression for /,, we start by setting Of /0t to zero in (15) and integ-ating twice, 
which gives: 





(20) a / . [ Siz yar dz 
z Das 0 
Expanding S as a series of orthogonal eigenfunctions gives: 
oo 
(20a) 52") = So an Fale") 
namQ 


and using orthogonality (18a) with (20a) gives: 


(20b) — i 32") 9 (2") dol”. 


Substituting (20a) into (20) gives 





20¢ 7 =f fy 7 (2) de" 
( ) Welz) 2 Duy 0 Xu An Jal )d daa 


Rearranging the order of summation and integration in (20c) gives: 





7 oo Xe 1 z' 
(20d) {.(2) = LS an | D / dale ide dz. 
aad z z's’ 40 
Now integrating (17a) twice gives: 


ze 1 x’ 
galz = / —y Xn Bee) dz di 
z ziz' 49 


or since Ap is a constant, 


(20e) dale). [ [ anl2t) ae" a2 
z 0 








(20f) Joo(2) = Doan Pat 
Substitution of (20b) into (20f) gives: 
(21) Jools) = So EL 3.2) a,(2) a? 


where S5(z2’) is positive (hence a source). 


Shulz and Lanzerotti (Ref. 45:39-40) point out that the directional flux 7(f;7r) = 
as f(y, P37) where f is the same distribution function as in this paper. Sh ilz (Ref. 47) 


points out that 7 must be integrated over all angles in the unit sphere in p-spz.ce to obtain 
the omnsdtrectronal flux J: 








cos @, 
(22 J snp? | f d(cos a). 
0 
Since 
sin’ a _ sin” Qo 
B Bo 
and 


it follows that 


and that j : 
—+(B/Bo) d 
d{cosa) = ee Enea 
V1 —y*(B/ Bo) 
Now since 


Bs 4 7) ee 
cos-a = 1 a =a Bw) 


if cosa = cosae, y"° = ye2; if cosa = 0, y* = Bo/B and (22) becomes 


t= any [” ~}(B/Bo) f-a(y?) 


Bo/B) \/1 — y2(B] Bo) 


which, upon combining terms and reversing the limits (-): 


(23) J= wf (B/ Bo) f d(y?) 
| y 


2 Vi—¥(B/Bo) 
Integrating (23) by parts 
[uae = w fod where u = f 


_ _(B/Bo) dy?) _ 


V1 — y?(B/ Bo) 


gives: 





and since 


will give (since the definite integral evaluates to zero, because f is 0 at y. by definition): 


ee , 
(24) Fm step | =) ) ots =). 


But, from (12), dz = —y T(y) dy. So 


dz _ -yT(y)dy _ -yT(y)dy _ —T(y) 
d(y2) dy?) ty dy 20 


Hence, (24) ultimately becomes: 


(25) opengl ae a7 .« 
J= —2xp? | Nhe (=) T(y) ae d(y). 
y a 


a 
€ 


Shulz then expands f as a weighted series of eigenfunctions: 


(26) fA ee inten tl gel=): 


rom 


This allows the flux J in (25) to be written: 


(27) (Bo/B) 
= 2xp° /- a 3 - ) Fela T(y) d(y*) 
(Bo/B) 
Sue L; af yi (4) T(y) F,(2) a} 


The following two expressions are derived in Appendix B. 


(28) 





(29) 


,_ V(2—¢)/2e (: - =) z =_— ( s yo 
In(2) Z JU (Kun) eee Ze Ju} Kun Ze 
‘ (1-o)/2 1~(@/2) (0/2) 
faa ! & Kun - o z 
+(2) 15(nee( 2) } a (1 =\(=) i 


_ The actual form of the flux equation awaits the choice of o, which also determines 
the dependence of D,, on z. Shulz uses o = O for purposes of comparisca of several 
computational methods. This choice is also attractive because it makes D,, independent 
of z (see 17(b)) which was a major assumption to begin with: 








(30) 1) 


This choice is made in the present study for the above reasons and also for the sake of 
simplicity in computations. 


lv 


With o = 0, the Bessel functios within eigenfunctions become order —1/, which can 
be simply written as: 


(31) J-1j2(R) =f 5 cos (R), 


(32) Ite) = -\/ = sin (J?) — a oan 


and the zeroes x, occur at (2n + 1)(x/2). 


This simplifies (32), since the second term will always be zero if J'_,/9(*n) is evaluated 
aka = (eme |)re). 


Equation (32) becomes 


£3 ] 
a 


ILialem= yl Tans aya sin (2n + 4) 


(33) J'_3(rn)= y/ Fin 1NR/2) (—1)"*?. 


The integral with respect to 2’ within (28) would thus be evaluated as: 





Bue » VV ot 
a : v f {2n+1)( wanti e/a 7=) N(z'/ ze) cos(( on +1) as 7) z 
: want ine/5) sin ((2n 1)7) = 








-/ V2/z_ cos ((2n + 1)(/2z¢) z’)dz! 


=sin ((2n + Sea 
((2n + 1)(7/2z2))(-1)"" 


(34) _ V2/2 sin((2n+ 1)(4/2)) 
((2n + 1)(x/2z-))(—1)""" 
The d, in (28) would be evaluated as (from (19)): 


ze 





(35) a = (2n aR v() Dee: 


we 


The flux equation (27), using (31), (33), (34), and (35), will then simplify to the following 
for o = 0: P 





(Bo/B) eel /zebstenlizee i (/2)) 
(35) J= ano f 
¥ 


Gguial ICES (- ve 
y? seesial 62 = 1)"( (x/2ze)* D 








: 2/2 


2 n+l wre 
EN (—1) Ze 
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- 


x(2n + 1)(x/2)(2/=.) ((2n : og 


& <e 


: (- aan inararay (e+ 5 2) 


= ——— os en Male ___) 
a((2n-+ I) /2)(z/2e))” yf samencTeNsTea 
n 9 a B ay 

Ene Utele) i= (4) T(y) ay?) | 


60 (Bo/B) 9 
2 ways IEE) a 


= 9 Ye? 





(x(2n + 1)(x/2)) (—1)"*" 


5] je ig ley se 
DzaN Ze m(2n + 1) 


panera (n+ OE =) 
Zz : 2 
uy vp (-sia((en+ ee n(2n + 1)(x/2)(2/ 20) 


7 cos ((2n + 1)(x/2)(z/ze)) 
a((2n + 1)(4/2)(z/ze))” /2/(x(2n + 1)(z/z)) 
; (2n + 1)(x/2) 


9 B Fd 
ss 2 jay, (=) T(y) a{y") i 
Equation (36) reduces to the following upon simplification and cancellation of terms 
i) (Bo/B) 
: on + 1)(x/2)(2/ze)) nee: 
Pee ony? / St 5(2/ze)sin ((2n + 1)(x/2)(2/ : eer ty) d 
a a Des,(2n + 1)*(x/2z-)*(—1)" Fn) AY) 
Bo/B - 
A sin((2n +1) r= )(en+ 0x] 


ft (2) \rw “o) 


c 





| 
(S> = 
a 
UU7y3 


The first term in Equation (37) represents a steady state value of the flux, and the second 
term represents a sum of eigenmodes. Bach A, (Ref. 45:162) has the form Gal Lea tr (EN 
As n increases, the r,’s become shorter; t.e., the higher eignmodes decay faster, so that 
eventually only the fundamental mode is left, which decays exponentially. At ¢ = oo, the 
fundamental mode has reached zero, and only the flux remaining is the steacy state flux, 
which corresponds to the natural environment background. 


The remaining substitutions in (37) are z., z, T(y), and y-: 


l 
VIB t= 18D) 


= Z(y)= (1 — y")T(0) - = (T(0) = T(1)](1 aa yon/ay) 
Z 


(33) fe 
Ze = Z(ye) = (1 — ye") T(0) - —{7(0) ~ r(yi( ~ ye) 


(as shown in Equation (14)) T(y) is given by Equation (11). 


The variable of integration in (37) is taken as y rather than y° for vimplicity of 
integration. 


The preceeding derivation has explicitly followed Shulz (1981). However, Shulz did 
not show how the differential-energy flux related to the integral-energy flu». Shulz and 
Lanzerotti (Ref. 45:163) state that the integral flux “will then scale as” the differential flux, 
without derivation. Since the AF WL trapped electron data base (Ref. 37) exists primarily as 
omnidirectional fuxes integrated over broad energy bandpass, the relation must be derived. 


It will be shown in the following and Lanzerotti (Ref. 45) are indeed -orrect: The 
differential-energy flux equation (37) will scale as the integral flux, and the energy dependence 
will be imbedded in the constants for the source and the cigenmodes. 


From Roederer (Ref. 39:86), the integral energy omnidirectional flux can be written: 


(39) Joe =f JdE 
E 


where J is the differential-energy flux, as in the previous equations, and Js<_ ir the integral 
energy flux above threshold energy E. 


A principal assumption will be that energy and pitch angle are both independent 
variables of the fux. Thus, (26) would really look like 


(39 a) f(z, £) = Fo, (E) Ja.{2) + AE 


nam () 


and hence in (25) Of /0z2 would look like 


(29 b } ar = Seo (E) Fooal2) + D> AnlE, Lit) a(2). 


rz () 


Using (27), (593) and (29b), Equation (39) becomes, at a given /-shell: 


2) 


(40) 
Jop=—trf pipes (EB) ap | Teo,(2) der dof p- 4B An,(B) | 7,(z) dz 
ZL ae & be | 


where Hi, (2) and 7,(z) represent all of the angular-dependent functions of z, 1’, or z shown 
in (25) through (27). 


Using the previously stated assumption (Ref. 45:162-163) that An, (E) are of the form 
an(Lye~*/'*-tE) 40) becomes: 


Joe=-2nf ye cB) 4B f }oo,(2) dz 


(4) ~ oxy. | page cme) an f jn(z) d 


rom () 


at a given /-shell. 


To “remove” the energy dependence of the decay-time terms, the expectation value 
of e~‘/(") is computed, assuming that the sccond term in (41) represents a listribution- 
function of the energies. 


(42) (<-¥ a) a, 


i p?a,(@) die 


Then (41) can be written, using (42): 
Joe =2n|- f- r*70,(6) 46) f ” Paq(e) ae 
(43) —2r | a én) e~t/(r0) ° g,,(z) dz. 
Yo [ fp rrentman] (om) fate 


ru =e Q) 1 


Now since it is assumed that the energy and angular shape functions are separable, 
the f,,,(#) must include only those parts of (28) which have energy ecceaa The only 
Lerms in (28) which could have any energy dependence are S and d, 


Shulz (Ref. 47:23) notes that, while- f(z) resembles go(z) in functional form, f,,(z) 
will coincide exactly with go(z) only if S$ is directly proportional to go(z). Fer simplicity, 
the assumption is made in this study that S is constant over the interval < z < 2z,. In 
physical terms, this really corresponds to a steady-state isotropic source. [t is known that 
the natural environment is not isotropic over cither Space or Lime, but the viriations are 
generally small compared to the variations from an injection by a nuclear evert. 


Now if S and \, are the only terms in (28) which may have energy depe idence, then 
, foo, (E) in (43) would correspond to 





=. S(Ie) 
De A 


rei) 7 


and since the only portion of the \, which may have any energy dependence {‘rom (19)) is 
D,:,, then f,,,(//) should be proportional to 


3(B) 
Deri) 


where D,,,(£) may be a constant in &. There may, of course, be some other functional 
energy-dependence included in f,,,(/) which is yet-undetermined. 


Equation (43) is generally of the same functional form as equation (37) if Lhe terms in 
square brackets in (43) correspond to the combined constant terms tn (37) as demonstrated 
above and as shown in Appendix D. Equation (43) also bears out Shulz and Lanzerotti’s (Ref. 
45:163) prediction that the integral flux scales with the differential flux; t.e., the “pitch-angle 
shape” dependence is invariant for differential- or integral-energy fluxes. 


Since it is impractical to evaluate an infinite number of terms in (37), only the terms 
up to n = 3 will be shown. Shulz and Lanzerotti (Ref. 45:163) state that the |igher modes 
vanish for n > 2. This in fact means that the higher modes decay so rapidly that modes 
higher that n = 3 should not be seen except at very early times. The final form of Equation 
(43), taken to n = 3, 18 shown in Appendix C. This is the form of the flux equ:.tion used tn 
the data analysis, except that fewer terms were used than are shown in Apperdix C. 


The general solution for the diflerential-energy flux (which, from (43), is directly propor- 
tional to the integral-energy flux), using Bessel functions of arbitrary order and with o- 
dependence explicitly shown, is given in Appendix D. 


Ya) 
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IV. ANALYSIS OF DATA 


-The AF WL Trapped Electron Data Base consists of a set of computer tapes of satellite 
counting data which have been assembled [from historical archives (such as the National 
Space Sciences Data Center, NSSDC) (Ref. 37). The data cover injections of elvctrons from 
nine high-altitude nuclear detonations shown in Table | (Ref. 6:6-2). 


The Data Base has been organized by Pfitzer into a coherent set of tapes which 
present, by satellite, net electron omnidirectional integral-energy (above threshold) fluxes as 
a function of time, 8, and L[-shell. The computed errors for each data point, in units of flux, 
are also given. The background which was subtracted to give net flux is also given (Ref. 
37:11, 188-189). 


Additional data have been collected by AF WL from plots published in early literature. 
These plots have been photographically enlarged, and the data points digitized and placed 
in computer files with formats similar to the Pfitzer data (Ref. 12:670-671, Ref. 31:646-468, 
Ref. 38:637-638). 


The satellites which provided the data, and their orbital parameters, are shown in Table 
Il (Ref. 44:41-45). 


The analysis of the data base required development of three computer >rograms by 
personnel at AF WL, in coordination with this author: 


(1) “Program DTABASE” to read raw data, organize it for processing, and output it 
into uniformly formatted files (Ref. 15). 


(2) A program to plot the raw data and to plot flux curves generated by the fitting 
program. 


(3) “Program Electrofit”, a fitting program to take data points from DTABASE at a 
given L-shell and energy and to fit them with the theoretical model developed in 
this study. It also generates flux curves using the fitted functions at predetermined 
times to compare with the raw data (Ref. 34). 


The first two programs were written with limited input from the author. The third was 
written under close and extensive coordination with this author, as discussed in this chapter. 


The DTABASE program, written by Miss Cherise Jarrett of Computer Sciences 
Corporation, reads the raw data files, in various formats created by Pfitzer (Ref. 37), satel- 
lite by satellite, and stores them in a three-dimensional array. The points are ‘hen ordered 
(within each energy group) by increasing time since burst. ‘To generate a file of points at a 
specified L-value, each data point L-value is compared to the desired Z and is recorded in 
the output file if a match occurs. Additional points are “created” by linear interpolation of 
two successive data points if they fall on both sides of the desired L and if the data points 
differ in time by no more than a specified value. The interpolation is perforried in L, 2, 
time, flux, and error. The output file can be restricted to a specified time “window” and to 
a specified range of B values, if desired; however, the data were analyzed over the full range 
of B and over all times which existed in the Data Base. The program has the additional 
capability to “B-average” data points at the same Lif their individual B values are within 
a user specified limit of each other and they are within the required time dillerence. This 
procedure has the effect of eliminating “double” points, and was used for the majority of 
this study with a limit of one percent. Another capability of the program is to interpolate 
to a specified B-value in a manner analogous to the /-interpolation described above; this 
capability was unused in this study (Ref. 15:1—2). 
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Name 


Alouette 1 


Alpha 
Upsilon 1 


Explorer 15 
Injun 1 
Traac 
Telstar 


Explorer 4 


TABLE 


Satellites and Their Orbital Parameters 


Launch 
Date 


9-28-62 


9-1-62 


10-27-62 


6-29-61 


11-15-61 


7-10-82 


7-26-58 


2-21-63 


10-86-58 


10-26-64 


10-23-59 


103.8 


105.6 


157.8 


110.1 


Perigee Apozee 


(mi.) = (mi.) 
620 633 
189 413 
194 10730 
534 631 
562 720 
593-3503 
163 «1372 


(Ref. 44:41- 45) 


In order to determine the optimum time difference to use for interpolation across the 
desired [-value, the DTABASE program was modified under direction of this author. The 
program now presents a summary table of the raw data points for each satellite, showing the 
frequency of occurance of time differences between successive data points from 0.1 minute 
to 5.0 minutes, by tenths of a minutes. The cumulative percentage of points with a specified 
time difference or less is also shown. These summaries are presented for the bursts and 
satellites examined as Tables Il] through XIII (Ref. 15:67-70). Perusal of these t::bles showed 
that the optimum time difference for interpolation was 2.2 minutes, primarily to include the 
maximum number of points in the Starfish Telstar file. This time difference wa: used for all 
data examined in this study. 


A final capability of DTABASE is to add data points to a specified L-value file by 
rounding points in L to the desired value. This procedure assumes that poin.s which are 
closer to the desired L than some specified limit may be considered to occur a! the desired 
L. Within the errors in satellite measurements, this procedure is reasonable, and in fact 
appears (by data comparison) to have been used by Roberts (Ref. 41:310), although he does 
not so explicitly state. This capability exists in the program, but was not used :n this study 
(Ref. 15). 


The plotting program was developed by Mr. John Burgio of AFWL to generate all the 
plota of raw data and fitted curves shown in this report. The program plots raw data points 
beginning at a specified start time since burst (in days) and covering the tine period in a 
specified time “window” (shown in Figure headings as “TW”). The L-value, equatorial B- 
value, and B-cutoff values are also shown. If the B-range of points is restricted, the minimum 
and maximim values of B are shown. The equatorial B value is computed from L by the 


lation 
relatio nage 


‘= 47 
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The B-cutoff value was initially computed by using the dipole formula (Ref. 39:55) 


1 
6 a = (Ba / Bee). 
; L3./4 — (3/1) eee 


However, it was later determined that this approximation might be less than optimum; so 
a table of L-values versus B.y~ was used and was linearly interpolated across the tabulated 
[-values to find the Beu, value. This table (Table XIV) was prepared by AF WL from the 
48-term Jensen-Cain model of the magnetosphere (Ref. 16). It was found that use of Table 
XIV produced Bey, (or yeut OF Zeut) Values more indicative of the global cutoff values. This 
table was used thereafter in both the plotting program and the fitting program (discussed 
below). This procedure is reasonable in view of the use of the same Jensen-C:1in model to 
produce the B-L coordinates for the Pfitzer Data Base (Ref. 37). 


The plotting program plots flux as a function of z = \/1 — (80/2). It also has the 
capability to plot flux data points versus time over all values of L and B, or flux versus 
energy at specified B or L. Next to the satellite name on each plot which has fitted curves 
is a number in equare brackets which corresponds to the serial number of the run of the 
fitting program which prodeuced the curves. 


This suthbor’s input to the plotting program was limited to supplying header information 
and specifying layout and format. 


The third program used in data analysis is “Program Electrofit,” written by Mr. Harry - 
Murphy of AFWL with the close coordination of this author (Ref. 34). Mlectro.it is a least- 
squares fitting program which uses the functions developed by this auther (.\ppendix C) 


g 
~ 


to fit data points which are output from DTABASE in order to compute the characteristic 
decay times (r,) and linear amplitude constants (a,). These decay times and constants are 
determined at the specific energy and L-value which are input with the data. 
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Table IJ 
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COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
ST ARF ISH 
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Table [V 
NO. DELTA 
IN BIN 
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(MIN) 


TIME BIN 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
ARGUS 1,2 


FOR ALL CHANNELS OF SATELLITE(S): 
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Table V 


(MIN) 


TIME BIN 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
RUSSIAN 1 


FOR ALL CHANNELS OF SATELLITE(S): 


FOR BURST: 
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Table VI 
EXPLORERI5S 
NO. DELTA TIMES CUMULAT IVE 
IN BIN PERCENT AGE 
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COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
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COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 


FOR ALL CHANNELS OF SATELLITE(S): 
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IN BIN PERCENT AGE 


TIME BIN 
(MIN) 


DODODDDCOCOCONNAAANANANMNDORRARARAR ARR DDDADAMADDADCVCVVVHH AAA ANO 


DOOSOCHSSSSSSSSSSOOS SS SSS SISSHIH ASIF GG ISFRPLEELLREK 


OO101O = O1SISIOIN = 016 91OIO'O MINIS 1111/0101010'010 00101010 00101010 7 eS aa 
0 al ¢ 


m 


HANNENARDOOANMNEINHARDFOAANTNHDRDOOKANTNHDARDAEOAANAM TNH OROGO 
DODDDDIDAD OH aan NANANANAANANANAMNAANMNNNeETeT TTT eT TO 
O 


ORANMENDRDOEOAANNENARDOOHRAMNENDRDOOHKANNITNDRDOOANNTNHOR GOGH 
CODD DDO D OO eae ae aH HAANNANANNANNAAOMMMMANANNNW TTT TTT TS 


Table VIII 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
FOR ALL CHANNELS OF SATELLITE(S): ALLOVETTE! 
FOR BURST: RUSSIAN 2 


TIME BIN NO. DELTA TIMES CUMULAT IVE 
(MIN) IN BIN PERCENT AGE 


As 1 
3107 78 
O 78 
409 89 
91 91 
7 91 
S2 92 
6) 92 
29 93 
25 94 
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i3 94 
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NO. DELTA TIMES 


Table IX 


(MIN) 


TIME BIN 


RUSSIAN 2 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
ALL CHANNELS OF SATELLITE(S): 


FOR BURST: 


FOR 
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TELSTAR 


IN BIN 


NO. DELTA TIMES 


Table X 


(MIN) 


TIME BIN 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
RUSSIAN 2 


FOR ALL-CHANNELS OF SATELLITE(S): 


FOR BURST. 
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CUMULAT IVE 
PERCENT AGE 


ALLOUETTE1 


IN BIN 


NO. DELTA TIMES 


Table XI 


(MIN) 


TIME BIN 


RUSSIAN 3 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
ALL CHANNELS OF SATELLITE(S): 


FOR BURST: 


FOR 
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Table XT] 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
FOR ALL CHANNELS OF SATELLITE(S): EXPLORER I5 
FOR BURST: RUSSIAN 3 


TIME BIN NO. DELTA TIMES CUMULAT IVE 
(MIN) IN BIN PERCENTAGE 

0Oo oO1 46 03 
0} 10-162 O 0 3 
0.2 03 77 08 
0.3 0.4 588 48 
0.4 O.5 1654 16 0 
0.5 0.6 2144 30.5 
0.6 0.7 1629 41 5 
0 7 70-68 970 48 0 
08 09 274 49.9 
09 1.0 194 Si. 2 
1-6. Tr 226 a2. 7 
Sly oie, 210 54.1 
1. 25s 293 56.1 
1 cee .4 273 5B 0 
14 1.5 474 61 2 
sy ler 442 b4 } 
{ieee 7 618 68 3 
17 1.8 313 70 4 
18 1.9 194 Ties 
1.9 2.0 131 72. 6 
2.Oye 1 82 Te 
2 ia 2 92 73 8 
22 2.3 74 74 3 
2.9 2.4 54 747 
24° 2&5 45 75.0 
2.5.) *2. 6 35 75 2 
. 2 cue 7 53 75.6 
27° 52.8 34 75.8 
ei 37 76 0 
29 30 40 76.3 
3) Ones 1 29 76 5 
3.4 eo 2 47 76 8 
32 33 52 77 a2 
393 64 b6 77 
34 35 53 78 0 
3.5 34 79 78. 5 
3.6 37 155 79 & 
a, aoe 91 80. 2 
3.8 3.9 127 81. 0 
3.9 F140 Vi a1 8 
40 4.1 @3 G2. 3 
ee a, 117 63 1 
4,2 9 43 87 83.7 
43 44 110 84.5 
44 4.5 73 85 0 
45 4.6 100 Q5. & 
4.6 4.7 76 BS 2 
4.7 4.8 63 Bb. & 
48 4.9 80 87 1 
49 5.0 59 Q7 5 
5. 0<= 1848 100.0 
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TELSTAR 
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Table XIII 


(MIN) 


TIME BIN 


COUNTS OF TIME DIFFERENCES BETWEEN NONZERO RAW DATA POINTS 
RUSSIAN 3 


FOR ALL CHANNELS OF SATELLITE(S): 


FOR BURST: 
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ro) 
OMANNENORDOOANMNEHGROCOANANEHORDOOAANTNHDORDOOANNTHORDONH 
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Table XTV 
Global Magnetic Cutoff Values 


[-Value Beutof sf 
1.0949 0.23901 
1.1152 22998 
1.1356 Bp Fo 
1.1559 22088 
babe aos 
1.2628 24017 
1.4817 .26770 
1.8331 roto 
2.3168 autos 
2.9330 33311 
3.6816 37001 
4.5626 40274 
5.9760 43178 

eos 15502 

8.0000 16951 


Note: The Table may be interpolated for other [-values (Ref. 16). 


The overall flow for the program is as follows (see Appendix I): 


(1) Use an initial guess for the fundamental decay time and the theoretical ratios 
of the higher-mode de«~y times to compute the starting decay tin es. The first 
guess at fundamental decay time is shown in Table XV as a function of L and 
energy. This table is interpolated in L and Ff for starting value. 


(2) Numerically compute the integrals of the functions shown in Appenclix C at cach 
data point using the initial guessed decay tiines. The integration 's performed 
by a 12-point Gauss/Legendre Quadrature at ench data point for ench function. 
Only the fundamental decay mode and two higher modes were cc mputed. 


(3) Compute the linear amplitude constants using a least-squares matrix-solving 
routine. Compute the standard error of the fit. The individual points are 
weighted by the reciprocal of their fractional errors, which have been normal- 
ized to an average value of one. 


(4) Optimize the initially computed yeue value by using a “Golden ‘Sicction Min- 
imim” function. This section attempts to let the data determin: the “best” 
cutoff value. If the data are not distributed down toward the less cone, the 
Golden Minimum function will provide an erroneous yey, value Hence, the 
function specifies that the computed y.y~ must be within +10 percent of the 
initial input value. If, in the iteration process, the function tries to exceed these 
limits, the initial guess of yeue is specified as the correct value. 


tcpeat the steps (2) and (3) with the optimum Yeue value. 


(6) To determine the best combination of decay times and amplituce constants, 
the decay times are changed by a “random-walk” method, either increased 
or decreased. The amount by which the decay times may randomly vary is 
determined by using 2 cumulative normal probability distribution about the 
initial value. 


The standard deviation used was estimated according to this author’s best guess of 
the accuracy of the tables and curves shown by Stassinopoulos (Ref. 53:31-32, 40-44) and 
West (Ref. 63:50-54). The Stassinopoulos tables do not always correlate with each other 
or with the Stassinopoulos curves, and the West curves do not always correlite with the 
Stassinopoulos in the /-regions of overlap. The initial fundamental decay times in Table XV 
were computed by this author by interpolation of the tables and curves in Reference 53 and 
by reading and extrapolation of the curves in Reference 63. The author’s estimates of their 
standard deviations are shown in Table XVI. Of particular note is the very low confidence 
in decay times at high L values (above 2.2) and at high energies (above 2.0 MeV). There 
have been few computations of decay times in these regions, and the results vary widely 
(Refs. 53, 63). 


The effect of using the cumulative normal probability distribution and the above stan- 
dard deviations is to restrict the “random-walk” of the decay times away from tlicir initial 
values. The farther away from the initial “mean value” that the random val ie generator 
places the new decay time, the greater the “push” that is created back toward -he mean for 
the next walk. The distance away, however, is defined in terms of the standa:d deviation, 
so that r valucs in which there is low confidence may vary by larger amounts than may r 
values in which there is greater confidence. The r’s nre also restricted by the co rstraint that 
To > Tt > Tt. 


With each change of decay times, steps (2) and (3) are repeated up to 96 times, in an 
attempt to get a “better” fit by reducing the standard error of the fit. 
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Table XVI 


Percent Standard Deviation of Initial 
Fundamental Decay Times (in Table XV) 


[-Range 1.02- L.21- L.4l- 1.81- = 
feo 1.40 1.80 2.20 ae 
Energy Range 
MeV 
E < 0.6 30 20 20 35 200 
0.8 << EF < 20 30 20 20 75 100 
ik > 3.0 20 20 50 100 100 


(estimated from Ref. 53:31-32, 40-44 and Ref. 62:50-54) 
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(7) If a better fit is found in less than 96 random walks, and if less than 32 “better” 
values have been found, steps (2), (3) and (6) are repeated up to 32 times, using 
the latest “best” decay times. 


(8) Up to five repeats of steps (4) and (5) may be carried out, but only if 32 “better” 
fits are found. Otherwise, only the first ycue optimization is performed. 


(9) The program will always end at the point where 96 random-walk3 have failed 
to reduce the error of the fit from the previous “best-fit”. The decay times, 
Yeur, 3nd amplitude constants are printed for each “better” fit, with the final 
values being those of the “best” fit. 


(10) After the best fit is found, the program uses the fitted constants and decay 
times to compute fluxes, at specified times, from the equator to the best cutolf 
value. These specified times correspond to default times plus 1/2 ‘ime window 
in {he program which plots raw data. Thus, the computed fluxer are used to 
generate curves which are plotted over raw data points at appro«imately the 
same times. 


Variations in the program were Lried in order to study the effects of the fitting process. 
The principal modifications were in the random decay-time variations. 


The first modification was to fix the ratios of the r’s to their theoretical values of 
1:1/9:1/25 (Ref. 48:16-19). The effect of this modification is to significantly reduce the 
degrees of freedom available to fit the data. 


The second modification was a compromise between the above and complete freedom to 
random walk. This variation allowed ro and r, to vary freely (subject to the usual contraint 
that m > 7), but r, and ro were fixed at their theoretical ratio of 25:9. The principal reason 
for this attempt was the very short time-span covered by some satellite/burst combinations. 
For example, the Russian 2 Burst covered only 41/2 days prior to the Russian 3 injection. 
If the data cover a time-span which is small compared to the fundamental decay constant, 
then the fitted value of that decay time may be suspect. Ilowever, if the time-s»an is not so 
small compared to the decay times of the higher modes, one should have more confidence 
in the fitted values of those decay time. Stated another way, the fit may not be sensitive to 
the fundamental, but may still be sensitive to the higher modes. 


In this analysis, not all L-values covered by the satellite in Table II were examined. 
Over 1000 plots of flux versus z were examined at various L-values for the satellite/burst 
combinations of Tables I and II. In many cases, the satellite orbits were cf such high 
inclination (toward polar) that data were not available close to equatorial B-values. Because 
of this, and because of the very large volume of data, those satellite/burst/L combinations 
which provided the most consistent coverage from the equator to the loss cone were chosen 
for study. 


Generally, the data fell into three regions of z-space (with some overlap). Region | was 
from the equator to the point where flux dropped off sharply into the loss cone (z = 0 to 
z = 0.6 or 0.8). Region If was the area where the flux “turned the corner” from a rather 
“horizontal” curve to arather “vertical” curve into the loss cone (z = 0.7 to 0.!'). Region TH 
was the loss cone region, where the flux levels dropped sharply to the cutoff vilue {z = .9 
for a typical L-value). 


Representative plots of raw data in the three regions are shown in Figures 4 through 8. 
These plots are typical in that data seldom cover all three regions in a short time-span. Gaps 
in the data appear because of satellite orbital coverage and because of lack of on-satcllite 
recorders. When a Satellite lacked a recorder, data were collected on!r as the satellite passed 
above the radio horizon of a ground tracking station. Figure 9 is an eight-hour plot of typical 
satellite coverage over all B-L space. 
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Figure 4(a). Raw Allouette Data Collected Only in Region MI (at burst tire). 
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Figure 4(b). Raw Alouette Data Collected Only in Region III (3 days post-burst). 
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Figure 5. Raw Telstar Data Collected in Regions I and Til (early time). 
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Figure 6. Raw Telstar Data Collected in Region LI (mid-range time). 
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Figure 7. Raw Telstar Data Collected in Regions I, Il, and I (mid-range ‘.ime). 
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Figure 8(a). Raw Telstar Data in Regions I and TI (late time). 
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Figure 8(b). Raw Telstar Data in Regions I, I], and Ii (late time). 
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Figure 9. Typical Raw Data Coverage for 8 Hours (B-L values not shown). 
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Figures 10 (a) through (g) shows raw data from Explorer XV over the time period of the 
Russian 2 and 3 bursts, at times corresponding to those shown in Figure 3, for comparison. 
It is apparent that fewer points are plotted from the AF WL data base in F gures 10 (a) 
through (g) than were plotted by Roberts in Figure 3 for the 1.9 Mc~ clectrons, even though 
the satellite and time periods of coverage were the same. These discrepancics are not resolved 
in this study, although the previously alluded-to L-rounding procedure nay account for the 
differences. 


The temporal progression of flux after one burst from one satellite (Telstar), at a typical 
L-value of 1.9, is shown in Figures 11 (a) through (h). Again, typical gaps in coverage of 
Regions I, H, and Ill are readily apparent over the 90 days shown. 


As previously noted, only fundamental and two higher eignemodes (fo, f1, f2) are used 
to fit the satellite data (Appendix C). The third eigenmode ( f3) is not used because Shulz and 
Lanzerotti (Ref. 45:163) suggest that the modes greater than two rapidly vanish. Initially, the 
steady-state function (f,,.) was used along with the fundamental, but this solution was found 
to be “competing” with the fit for the fundamental, in the sense that the linear constants 
tended to be roughly equal in magnitude and tended to alternate in sign. Shulz (Ref. 47:22- 
23) indicates that the f,, and fp functions should resemble each other ia shape. Also, Pfitzer 
(Ref.37) has subtracted the background from all flux data, and the background should 
approximate the steady-state, f.., solution. For these reasons, the steady-state solution was 
not used in the data fit. 
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Figure 10(a). Raw Explorer 15 Data at Times Corresponding to Figure 3°. 
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Figure 10(b). Raw Explorer 15 Data at Times Corresponding to Figure 3. 
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Figure 10{c). Raw Explorer 15 Data at Times Corresponding to Figure S. 
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Figure 10fd). Raw Explorer 15 Data at Times Corresponding to Figure 3. 
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Figure 10(e). Raw Explorer 15 Data at Times Corresponding to Figure 3. 


60 


INTEGRAL ELECTRON FLUX (FLECTRONS/CM*)/SEC 


BURST : RUSSIAN 2 


SAT : EXPLORERI5/EXPLORERI5 


TIME : 47 (dy) 19 (hr) 22.06 (min) 47.807(days) 
Pec Oee beqe..04908CUT > 0.479 TW (days): 2.0 
; Bmin (Gauss): 0.001 Bmax (Gauss) : 1.000 
10 (-CEGEND 
O = .500 MEV 

10 

10° t 

10° } 

10° 

10° 

10° Ses alee 

0 


05 
X = (1-Beq/B) / 


Figure 10(f). Raw Explorer 15 Data at Times Corresponding to Figure 3. 
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Figure 10(g). Raw Explorer 15 Data at Times Corresponding to Figure 3. 
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Figure 11(a). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(b). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(c). Temporal Progression for Raw Telstar Data from | to 90 Days Post-Burst. 
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Figare 11(d). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(e). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(f). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(g). Temporal Progression for Raw Telstar Data from 1 to 90 Days Post-Burst. 
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Figure 11(h). Temporal Progression for Raw Telstar Data from 1 to 90 Day: Post-Burst. 
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Panel) sb) DISCUSSION 


Figures 12 and 13 show the typical temporal progression of raw data at two representative 
L values, L = 2.3 and L = 2.4. Included on those figures are the flux curves calculated from 
the amplitude coefficients and decay times fitted to that data by Electrofit. As previously 
noted, the curves shown on any given figure may vary slightly in time from the individual 
data points, since the curves represent only the one time-instant corresponding to the start 
time plus one-half the time window of the plot. Figures 12 and J3 represent maximum 
freedom of the fit, since each r is free to vary independently of the other r’s. 


Figure 14 shows the same temporal data progression with different fitted curves. The curves 
in this figure were calculated under the restriction that the ratios of the r’s were fixed at 
the theoretical values of 1:1/9:1/25; hence, only ™ was really varying in a random sence. 


Figure 15 shows similar data progression and fitted curves where the r’s were allowed to 
freely vary, subject to the constraint that m > ™ > mm. Also, this was an initial run with 500 
tries for convergence, vice 96 tries. [n addition, the computed flux was used in the iteration 
process to eliminate non-physical (negative) solutions. The solution was restricted to being 
positive at four different times and at five different equatorial pitch angles. Altiough this is 
an initial run with only a small amount of data used, it indicates that the reconimendations 
in Chapter VI should be pursued. 


Table XVII shows the fitted linear amplitude coefficients and exponential decay times (in 
days) at constant Energy and /[-shell for a representative subset of the data examined in 
this study. The sheer volume of resultant fits and calculated curves made it impossible to 
include all results in this report. The valnes shown in Table XVII were computed with the 
r’s allowed to vary freely. , 


Figures 12 and 13 are considered representative of the results of the fitting model. It is 
apparent from examination of Figures 12 and 13 that the pitch-angle dillusion theory 
and flux computation method developed in this report are consistent with experimental 
measurements. Higher eigenmodes are readily discernible at the lowest energy in the Figures 
or the first five days, and are not visible at later times. These higher eigeumodes are not 
discernible in the plots at the higher energies; Table XVII shows that these modes are 
present, but at much lower levels, in some cases. The calculated curves show 1 reasonable 
fit to the data even at late times, consistent with the assumption of exponential decay. 


Examination of Figures 14 and 15 shows that some inconsistencies remain to be resolved in 
the use of this fitting method. The caluclated flux curves for the low-energy data points also 
appear to be slightly lower than optimum, while the high-energy data appears to be fitted 
better by the model. This may be partially explained by the existence of many lower data 
points in the vicinity of the loss cone which tend to “drag down” the overall fit. Any change 
in the cutoff valuc significautiy affects the Gt by shifting the “vertical” portion of the curve 
toward or away from those points. However, this cannot account for the entire difference 
of the curves from the data. The energy dependence may be treated less thin optimally 
in the choice of functional solution made in this study. If the chosen function solution is 
not the correct one, it may be nevertheless close enough to correct to give a r2asonable fit 
at high energy, but be significantly in error at lower energy. A different functional solution 
mav fit the lower-:nergy groups better and yet still lit the high-energy elect-ons as well. 
Recommendation for the use of other functions is made in Chapter VI. 
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Figure 12(a). Temporal Progression of Raw Data and Fitted Flux Curves at b = AST): 


Russian 3 Burst (r, freely varying). 
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Figure 12(b). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 


Russian 3 Burst (7, freely varying). 
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Figure 12(c). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (r, freely varying). 
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Figure 12(d). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (ty freely varying). 
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Figure 12(e}. Temporal Progression of Raw Data and Fitted Flux Curves at 7, = 2.3 for 
Russian 3 Burst (7, freely varying). 


10 


10 


INTEGRAL ELECTRON FLUX (ELECTRONS/CM*)/SEC 


Be Rote RUSSIAN Seen LORE IS 1045) 
TIME: 4 (dy) O(hr) 0.00 (min) 4.000(days) 
be 2.50) Bea: -025 BCUT: 0.506 TW (days): 1.0 


10 


10 





Bmin (Gauss): 0.001 Bmax (Gauss) : 1.000 


10 


LEGEND 10 
O = 500 MEV 
O = .500 MEV 
& = 1.90 MEV 
3 “< = 1.90 MEV ; 
HOS = 2.90 MEV 10 
VY = 2.90 MEV 
10 10 
10° 
10° 
10 
10. te a AY ce 10. 


OF) 
X= (1-Beq/8)'”” 
Figure 12(f). Temporal Progression of Raw Data and Fitted Flax Curves at /, == 2.9 for 
Russian 3 Burst (rq freely varying). 
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Figure 12(g). Temporal [rogression of Raw Data and Fitted Flax Curves at I, = 2.3 for 


jLussian 3 Burst (rp freely varying). 
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Figure 12(h). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (rq freely varying). 


10 


INTEGRAL ELECTRON FLUX (ELECTRONS/CM“)/SEC 


BURST : RUSSIAN 3 SAT : EXPLORER15 [045] 
TIME : 60 (dy) O (hr) 0.00 (min) 60.000(days) 
L:2.30  Beq:.025 BCU: O S05 IW (das: S50 

Bmin (Gauss) : 0.001 8max (Gauss) : 1.000 


10 


10 





LEGEND 
O = .500 MEV 
O = 500 MEV 
10 
10° 
10 
10° | | 
ie) 
10 
0 


05 an 
X = (1-Beq/B) / 


Figure 12(i}. Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (rp, frecly varyiog). 
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Figure 13(a). Temporal Progression of Raw Data and Fitted Flax Curves at 1 = 2.4 for 
Russian 3 Burst (r, freely varying). 
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Figure 13(b). Temperal Progression of Raw Data and Fitted Flux Curves at L = 2.4 for 
Russian 3 Burst (r, freely varying). 
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Figure 13(c). Temporal Progression of Raw Data and Fitted Flux Curves at 1 = 2.4 for 
Russian 3 Burst (7, freely varying). 
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Figure 13(d). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.4 for 
Russian 3 Burst (rT, freely varying). 
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Figure 13(e). Temporal Progression of Raw Data and Fitted Flux Curves at = 2.4 for 
Russian 3 Burst (7, freely varying). 
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Figure 13(f). Temporal Progression of Raw Data and Fitted Flax Curves at 1 = 2.4 for 
Russian 3 Burst (tr, freely varying). 
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Figure 12(g). Temporal Progression of Raw Dats and Fitted Flux Curves at i = 2.4 for 
Russian 3 Barst (7, freely varying). 
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Figure 13(h). Temporal Progression of Raw Data and Fitted Flux Curves at 1 = 2.4 for 
Russian 3 Burst (rt, freely varying). 


gg, 


INTEGRAL ELECTRON FLUX (ELECTRONS/CM*)/SEC 


BURST : RUSSIAN 3 SAT : EXPLORERI5 [052] 
TIME: O(dy) O(hr) 0.00 (min) 0.000(days) 
ic COmmpeo | O25.,8CUT 70.506 TWid@ays): 0.5 
Bmin (Gauss) : 0.001 Bmax (Gauss) : 1.000 


10 


LEGEND 


gaoxXK bOOD 


— 
= 
— 
=> 
— 
— 
== 
—— 
= 
= 
= 
= 





a 
X = (i-Beq/B) fe 


Figure 14(a). Temporal Progression of Raw Data and Fitted Flux Curves at LD = 2.3 for 
Russian 3 Burst (r, fixed 1:1/9:1/25). 
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Figure 14(b). Temporal Progression of Raw Data and Fitted Flux Curves at 1 = 2.3 for 
Russian 3 Burst (rp fixed 1:1/9:1/25). 
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Figure 14(c). Temporal Progression of Raw Data and Fitted Flux Curves at & = 2.3 for 
Russian 3 Burst (tp, fixed 1:1/9:1/25). 
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Figure 14(d). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (7, fixed 1:1/9:1/25). 
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Figure 14(e). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (r, fixed 1:1/9:1/25). 
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Figure 14(f). Temporal Progression of Raw Data and Fitted Flux Girvessate — 2-omler 
Russian 3 Burst (rp, fixed {:1/9:1/25). 
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Figure 14(g). Temporal Pregression of Raw Data and Fitted Flux Curves at L = 2.3 for 
Russian 3 Burst (rp fixed 1:1/9:1/25). 
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Figure 15{a). Temporal Progression of Raw Data and Fitted Flax Curves at 1, == 2.4 for 
Russian 3 Burst (rt, freely varying with physical contraints). 
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Figare 15(b). Temporal Progression of Raw Data and Fitted Flux Curves at 2 = 2.4 ‘or 
Russian 3 Burst (rz freely varying with physical contraints). 
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Figure 15(c). Temporal Progression of Raw Data and Fitted Flux Curves at 1 = 2.4 for 
Russian 3 Burst (r_ freely varying with physical contraints). 
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Figure 15(d). Temporal Progression of Raw Data and Fitted Flux Curves at 2 = 2.4 for 
Russinn 3 Burst (rz freciy varying with physical contraints). 
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Figure 15(e}. Temporal Progression of Raw Data and Pitted Flux Carves at 1 = 2.4 for 
Russian 3 Burst (rq freely varying with physical coatraints). 
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Figare 15(f). Temporal Progression cf Raw Data and Fitted Flux Curves at L = 2.4 for 
Russian 3 Burst (7, freely varying with physical cootraints). 
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Figure 15(¢). Temporal Progression of Raw Data and Pitted Plax Curves at L = 2.4 for 
Russian 3 Burst (7, freely varying with physical contraiuts). 
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Figure 15(h). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.4 for 
Russinn 2 Burst (7, freely varying with physical contraints). 
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Figure 15(i). Terporal Progression of Raw Data and Fitted Flux Curves at L = 2.4 for 
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Figure 15(j). Temporal Progression of Raw Data and Fitted Flux Curves at L = 2.4 for 
Russian 3 Burst (ra freely varying wth physical contraints). 
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Figure 15(k). Temporal Progression of Raw Data and Fitted Flux Curves at , == 2.4 for 
Russian 3 Burst (7, freely varying with physical contraints). 
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Figure 15(1). Temporal Progressien of Raw Data and Fitted Flux Curves at /, = 2.4 for 
Russian 3 Burst (7, freely varying with physical contraints). 
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Figure 15(m). Temporal Progression of Raw Data and Fitted Flax Curves at 1 = 2.4 for 
Russian 3 Burst (rq freely varying with physical contraints). 
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The short time-spans covered by most of the data sets in the Trapped Electrcn Data Base 
have not allowed adequate determination of m. The values of m presented in Table XVII 
are suspect for this reason. However, more confidence is held in the values of 7, and 19 
computed by this study. In most data sets studicd, the time period of the data was long 
enough compared to the values of 7 and 7; to compute reasonable values. More consistency 
between data sets was found for 7, and ro than was found for ro. Of course, then tm may be 
simply computed by multiplying 7, by the theoretical ratio corresponding to the solution. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 


The pitch-angle diffusion theory and flux calculation method developed in this study 
are reasonably consistent with experimental data in most cases. The method si ows promise 
as a way of calculating improved decay times over a broad region of the mazsnetosphere, 
particularly for higher energy clectrons. Additionally, a comprehensive summary table of 
the best available literature values for ro has been prepared (Table XV). A complete table 
covering such a broad range of energy and L has not heretofore been compiled. 


Time and computer limitations have precluded examination of the entire data base. 
However, the overall efficiency of the model developed in this study has been demonstrated. It 
is expected that improved predictions of fluxes and improved inputs to the AFW', SPECTER 
Codes (Ref. 7-9) will be available as a result of this work. This will enhance tie capability 
for calculation of satellite operational environments, and will improve the abili.y to predict 
satellite survivability and vulnerability. 


Recommendations 


The following recommendations for further study are proposed to improve the model 
developed in this study, and to resolve the problem of non-physical! solutions for some data 
sets. At the direction of this author, several of these recommendations are presently being 
implemented at the Air Force Weapons Laboratory. 


(1) Incorporate an estimation of the error in the magnetic field value, B, of 
individual data points. This corresponds to a “horizontal” error bar on the 
plots such as Figure 12. It is known that errors exist in the computation of the 
B value associated with each data point, and this error would prove significant 
in the region of the loss cone. Improved weighting of the data points with the 
combined errors should improve the fit. 


(2) Utilize the L-rounding procedure of Chapter ITV to add data points to each L- 
value considered. This should aid in filling gaps in the data, whic necessarily 
will improve the fit. 


(3) Discard the highest eigenmode if the fitted value of its decay tim» is less that 
0.25 day, and fit the data with only the fundamental and one higher mode. 
A value of 0.25 day is not reasonably fitted within the errors in the Trapped 
Electron Data Base. This mode would be essentially lost within th» first day in 
any case. 


(4) Perform an “energy-scaling” of the data to combine data from different satel- 
lites. This involves assuming some form of energy spectrum, such as e~!4//e) | 
including an estimate of the error in the assumed value of Ko. The new flux 
would be simply the old flux multiplied by a factor e~(Pree—Faral/Ee. If the 
error in the old flux were 50 percent (a value consistent with the Data Base) 
and the error in Ho were 25 percent, and if Eo = 1.1 (a fission spectrum), then 
the resultant error in the shifted flux would be less than 1 percent different 
from that of the old flux, when shifting from 660 KeV to 500 KeV. Thus, it 
should be possible to combine data sets and fill the gaps in data to achieve a 
better fit. However, the limitation of the assumed spectrum still exists. 
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(5) Utilize only late-time data where it is available, and fit only the fundamental 
mode to it, under the assumption that higher modes will have decayed away. 
Then, with the fundamental mode constrained, go back to early time data and 
attempt to fit higher modes, if any. This recommendation is of limited value 
because only two data sets cover a significant period of time: Starfish/Telstar 
and Russian 3/Ixplorer 5. However, it may provide a better value of rp than is 
presently available at some JL values. 


(6) Revise the fitting functions to other Bessel function solutions than the Joise 
solution and compare results. This was performed only on one other functional 
solution, Jo, with only one small data set, and the results showed a0 improve- 
ment over the present study. However, to be certain that the auhor has not 
simply chosen the wrong functional form, more functional forms must be in- 
vestigated. 


(7) In the fitting program, constrain the directional flux to physical reality at each 
iteration (i.e., compute flux at each iteration and reject non-physical solutions 
at each iteration, rather than after converging on a final solution). This should 
converge to a physically real solution if it converges at all; however, significant 
increases 1D computer processing time may accrue. 


(8) Increase the number of tries for a solution from 96 to 200 and reduce the number 
of outer loops from 32 to 10 or 15. This may improve the convergence of the 
least squares error, if non-physical solutions are also rejected at cach iteration. 


(9) Compare the background subtracted from the flux with the latest available 
natural environment, to ascertain the correctness of the background subtrac- 
tion. If significant differences are found, the background should be re-added to 
the fiux and the natural environment subtracted. 


(10) A suggestion proposed by Professor D.G. Shankland, Air Force Institute of 
Technology, involves a significant modification to the fitting method. This 
method involves an interpolating function, or measure of “roughness” of the 
fit, which may be made arbitrarily smooth. The advantage of this method is 
that no limiting assumptions need be made about the data. The chief disad- 
vantage is that the method is a strictly mathematical method, with no physical 
constraints. Although the method would involve extensive reprogramming, it 
promises worthwhile information about the Mux curves, and shonld be tried. 
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APPENDIX A 


Demonstration of Zero Diffusion Current 
at the Magnetic Equator for Half-Integer 
Bessel Function Solution 


To demonstrate zero diffusion current, onc must show that 
lim,—o De 7',(2) 1 


for the o = O case. We Know that from (31): 


Juj2(2) = V2/(r2) 008 (2) 


and from (16): 
De — Digs: 


and from (18): 


V 2/2 2 F Az 
oe n+ 1)(x/2(z/z) ((2n ea =) 


In(z) = JI" en) ) 


= i slices abate erase) =m cos ((en +1) 2) 


J (Ra) 


Hence, the derivative of the cigenfunction: 


~| Seven + ra sin((2n +, eT 


oR) ee 


Finally, the quantity which must go to zero is: 


Pes V2] 20 J 1/ (2 (2n + 1) ((2n + 1) “tea sin((on bis a 


Dea. = 
In J"(kn) 


The term in square brackets is a finite constant for any n, and clearly the sine goes to zero 


as z goes to zero. Hence, the limit is zero. 
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APPENDIX B 


Derivation of Figenfunction Derivatives f,, and 7’, 


To find the derivative of the steady state function, f.,, one starts with (21): 


Too = ode fe) ale a2) 


If the 4, are assumed (from (19)) independent of z, then the series can simply be 
differentiated term-by-term, since the definite integral is simply a constant for each value 


of n. | 
H teen £| 53> gale) [5 
fo.{2)= aps Dies i 5 In(=) ae) 


') d(z’). 


tJ2 
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To find the derivative of the eigenfunction g,, one starts with (18): 


vie= eM : i (- ayo") 





Gn(z) = ae es 


SS 


- 
“¢ 


Using the product rule for differentiation: 


o _ ~ \a~lett)/2 ~ \(l=(¢/2)) 
Bp = ey tye (2) 
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APPENDIX C 


The Specific Solution of the Integral-Energy 
Omnidirectional Flux Equation for the BES B32) 
Bessel Function Solution 
The integral-energy omnidirectional flux Mquation (43), which is the energv-integrated 


form of Equation (37), is represented as the steady stale term, f.., plus the sum of exponen- 
tially decaying eigenmodes fg through fs: 


J>e(y = VBo/B) = foo + fo + hi + fe + fs 


In this appendix, each term is shown, as derived from Equations (37) and (43). Note 
that the eigenmode sum should be infinite, but only the fundamental mode and three higher 
modes are shown. 


The constants in each term are defined as follows (Ref. 39:55): 
ye = f1/( V4— (3/0) 
T) = 1+ 1G + vil) (evs)|  1,3801730 


— (xv2)/6 = 0.7404805 
K = (4/11)[T> — Ty} = 0.2328155 
7, = The K(1 = yi) 


The “fitting constants” are defined as in Equation (43) or in Appendix D: 
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These fitting constants are really fuctions of energy and /-shell. 
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The ratios of ro, 7,, T2, and 73; should be (as computed from the cigen value ratio of the 
zeroes of the function): 


To oz Or; aa (1/9) r0 
To Fe 25 ro or alternatively Fo (1) oir 
To FS AIrs r; =(1/49)% 


APPENDIX D 


The General Solution of the Integral-Energy 
Omnidirectional Flux Equation for Bessel 'unctions 
of Arbitrary Order, as a Function of Sigma 


Substitution of (18), (28), and (29) into (27) completes the general expression for the 
differential energy omnidirectional flux, J, with Bessel functions of order 
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where o is some number less than two: 
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The analogy between the above equation and Equation (43) is seen if: 


| | D eens) dE = / o°S/Diz. di 
E L Ee 
and if: ah . a 

"a n(F SOG =| a a mea 
I. Qa i (e— €/(Tn)}p ae Les 


x 
== [ yw an(Hle7 fl EN ap 


vw ft 


Henec, (43) is simply the energy-intesrated form of (37). 


APPENDIX 


Flow Chart for Program Electrofit 


Initial guess t) (and 71, ™). 






Numerically integrate functions at 
each data point. 









Least squares fit for the constants: 
ve. 


_——— oOo lc —_— ee ii ->ll SSS | Optimize Yerut 


Compare fit against previous “best” 
fil using least squares residual. better 


Update “best” r’s and a's. 


no better 


Check for 96 tries with no im- 
provement 
Check for 32 better fits | yes 


no 








Random walk r's 





Print results. 


Compute fluxes. 
Plot computed fluxes with raw 


data. 
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